
%%%%%%%%%%%%%%%%%%%
%%% Table 7 %%%
%%%%%%%%%%%%%%%%%%%

clc
clear all

% variable set-up

params.alpha=0.5;
params.eta=6;
params.r=0.02;
params.pi=0.001;


%%%%%%%%%%%%%%%%%%%
%%% External Innovation Cost 10% higher
%%%%%%%%%%%%%%%%%%%

%%% Counterfactual f=0

    load calibration_outcome_f0.mat

    params.lambdaE=xx_cal(1);
    params.lambdaS=xx_cal(2);
    params.lambdaI=xx_cal(3);
    params.lambdaC=xx_cal(4);
    params.f=0;
    params.xiE=xx_cal(5)*1.1;
    params.xiI=xx_cal(6);
    params.lambdae=params.lambdaE;
    params.xie=xx_cal(7);

    guess_sim=0.1*ones(10,1);
    options = optimset('MaxFunEvals',10^5,'MaxIter',10^5,'TolFun',1e-30);
    [xx_sim,fval_sim] = fminsearchbnd(@(x) equation_difference(x,params),guess_sim,[0 0 0 0 0 0 0 0 0 0],[],options);
    
    xx_eqm_counterfactual=[xx_sim(2), xx_sim(5), xx_sim(8), xx_sim(7)];

    percent_difference=(xx_eqm_counterfactual-xx_eqm)./xx_eqm*100
    

%%%%%%%%%%%%%%%%%%%
%%% Baseline f=17.6
%%%%%%%%%%%%%%%%%%%

    load calibration_outcome.mat

    params.lambdaE=xx_cal(1);
    params.lambdaS=xx_cal(2);
    params.lambdaI=xx_cal(3);
    params.lambdaC=xx_cal(4);
    params.f=xx_cal(5);
    params.xiE=xx_cal(6)*1.1;
    params.xiI=xx_cal(7);
    params.lambdae=params.lambdaE;
    params.xie=xx_cal(8);

    guess_sim=0.1*ones(10,1);
    options = optimset('MaxFunEvals',10^5,'MaxIter',10^5,'TolFun',1e-30);
    [xx_sim,fval_sim] = fminsearchbnd(@(x) equation_difference(x,params),guess_sim,[0 0 0 0 0 0 0 0 0 0],[],options);
    
    xx_eqm_counterfactual=[xx_sim(2), xx_sim(5), xx_sim(8), xx_sim(7)];

    percent_difference=(xx_eqm_counterfactual-xx_eqm)./xx_eqm*100
    






%%%%%%%%%%%%%%%%%%%
%%% All Innovation Cost 10% higher
%%%%%%%%%%%%%%%%%%%

%%% Counterfactual f=0

    load calibration_outcome_f0.mat

    params.lambdaE=xx_cal(1);
    params.lambdaS=xx_cal(2);
    params.lambdaI=xx_cal(3);
    params.lambdaC=xx_cal(4);
    params.f=0;
    params.xiE=xx_cal(5)*1.10001;
    params.xiI=xx_cal(6)*1.10001;
    params.lambdae=params.lambdaE;
    params.xie=xx_cal(7)*1.10001;

    guess_sim=0.1*ones(10,1);
    options = optimset('MaxFunEvals',10^5,'MaxIter',10^5,'TolFun',1e-30);
    [xx_sim,fval_sim] = fminsearchbnd(@(x) equation_difference(x,params),guess_sim,[0 0 0 0 0 0 0 0 0 0],[],options);
    
    xx_eqm_counterfactual=[xx_sim(2), xx_sim(5), xx_sim(8), xx_sim(7)];

    percent_difference=(xx_eqm_counterfactual-xx_eqm)./xx_eqm*100
    

%%%%%%%%%%%%%%%%%%%
%%% Baseline f=17.6
%%%%%%%%%%%%%%%%%%%

    load calibration_outcome.mat

    params.lambdaE=xx_cal(1);
    params.lambdaS=xx_cal(2);
    params.lambdaI=xx_cal(3);
    params.lambdaC=xx_cal(4);
    params.f=xx_cal(5);
    params.xiE=xx_cal(6)*1.10001;
    params.xiI=xx_cal(7)*1.10001;
    params.lambdae=params.lambdaE;
    params.xie=xx_cal(8)*1.10001;

    guess_sim=0.1*ones(10,1);
    options = optimset('MaxFunEvals',10^5,'MaxIter',10^5,'TolFun',1e-30);
    [xx_sim,fval_sim] = fminsearchbnd(@(x) equation_difference(x,params),guess_sim,[0 0 0 0 0 0 0 0 0 0],[],options);
    
    xx_eqm_counterfactual=[xx_sim(2), xx_sim(5), xx_sim(8), xx_sim(7)];

    percent_difference=(xx_eqm_counterfactual-xx_eqm)./xx_eqm*100
    
